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Although quantum computers have the potential to efficiently solve certain problems considered 
difficult by known classical approaches, the design of a quantum circuit remains computationally 
difficult. It is known that the optimal gate design problem is equivalent to the solution of an 
associated optimal control problem, the solution to which is also computationally intensive. Hence, 
in this article, we introduce the application of a class of numerical methods (termed the max-plus 
curse of dimensionality free techniques) that determine the optimal control thereby synthesizing 
the desired unitary gate. The application of this technique to quantum systems has a growth in 
complexity that depends on the cardinality of the control set approximation rather than the much 
larger growth with respect to spatial dimensions in approaches based on gridding of the space, used 
in previous literature. This technique is demonstrated by obtaining an approximate solution for the 
gate synthesis on SU(4)- a problem that is computationally intractable by grid based approaches. 
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I. INTRODUCTION 

The advent of Shor's algorithm [T] demonstrated the 
potential for processors based on quantum operations to 
perform certain computational tasks exponentially faster 
than those limited to to using classical operations. There 
has been much work devoted to solving the following 
problem of special interest: to determine the bounds on 
the number of one and two qubit gates required to per- 
form a desired unitary operation-termed the gate com- 
plexity of the unitary. Yet, the explicit design of quantum 
algorithms has remained a challenging task. 

One approach to this task of constructing an optimal 
circuit was highlighted in [5] where it was shown to be 
equivalent to finding a least path-length trajectory on a 
Riemannian manifold. This insight opened up the study 
of quantum circuit complexity to the use of tools from 
optimal control theory. 

In [3] the method of dynamic programming was in- 
troduced to solve the control problem associated with 
quantum circuit complexity. The numerical computa- 
tions of solutions using this technique proceeded via a 
widely used grid (mesh) based iteration approach 
that requires the generation of a mesh in the region of 
the state space over which the solution is sought. This 
approach however, leads to the following issue. A grid 
(assumed for simplicity to be regular and rectangular) 
with K points along each of the N dimensions has K N 
grid points, over which the solution must be propagated 
during each iteration. In addition, the dimension of an 
n qubit quantum system grows exponentially (as 4" — 1) 
thereby leading to a similar exponential growth in mem- 



ory and time requirements. This large growth in the re- 
sources required, arising from growth in the dimensions of 
the system, is termed the curse of dimensionality (COD). 
It renders the direct application of mesh based solution 
techniques unfeasible for systems larger than SU(2) due 
to the large memory (in terabytes) and time (in centuries) 
required to solve problems of these dimensions via such 
methods. 



In |j a COD-free technique was introduced for 
problems in Euclidean space. In this article we adapt 
these methods for quantum systems. Due to the struc- 
ture of the control problem that we consider, we do 
not completely eliminate the COD. However we have a 
much more manageable growth related to the number 
of elements in the discretized control set used. This is 



managed via a pruning approach described in Sec. IV 
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The computational time of the resulting algorithm grows 
much slower than that in mesh based methods, thereby 
bringing us closer to the numerical study of larger sys- 
tems. One particular application of interest for the 
numerical methods developed is the determination of 
whether a given unitary U in an n-qubit system can be 
approximately synthesized in an efficient manner (with 
respect to the growth in n) in a given time T. 



The paper is structured as follows. Sec. |TT] gives a 
brief introduction to the relevant concepts in quantum 
complexity and optimal control. We then introduce the 
reduced complexity algorithm in Sec. |III[ and in Sec. IV 
highlight the complexity growth in the application of this 
method and its management. The algorithm is then ap- 
plied in Sec. [V] to the two qubit optimal gate synthesis 
problem on SU(4) - a problem in 15 dimensional space. 
In Sec. |VI| we conclude with comments on various aspects 
of the technique introduced in this article. 
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II. PRELIMINARY CONCEPTS 

In this section we recall the notion of gate complexity 
and introduce the cost function for an associated control 
problem as in [31 |S] . 

A. Gate complexity and control 

In quantum computing an algorithm operating on an 
n qubit system can be represented as an element of the 
Lie group SU(2 n ) (denoted in this article by G) and 
is termed a unitary. Every such unitary can be con- 
structed by a sequence of available elementary unitaries 
Ux, U2 • ■ ■ U n . In practice, we synthesize a unitary Uq that 
approximates a desired computation Uq with a required 
accuracy e (i.e. ||Z7o — Uq\\ < e, where || • || denotes the 
standard matrix norm). This leads to the notion of ap- 
proximate gate complexity G{Uq, e) which is the minimal 
number of one and two qubit gates required to synthesize 
Uq up to an accuracy of e without ancilla qubits [S]- 

Related to the gate synthesis problem is an optimal 
control problem (described below) on G, such that the 
approximate gate complexity scales equivalently up to a 
polynomial in the optimal cost function for the control 
problem. This equivalence motivates the solution of the 
associated control problem. 

We now describe the control problem and recall the 
solution process, via the dynamic programming principle, 
used in [3j. 

B. System Description 

The system dynamics for the gate design problem is 
given by: 

^=-i{Y.Vk(t)Hk}U, UeG (1) 

fe=i 

with control v (such that v(t) G R M , Vi > 0) and 
an initial condition U(0) — Uq. For the class of prob- 
lems considered, v is taken to be an element of the set 
of piecewise continuous functions having a norm bound 
ll u (')ll = 1 (where || • || denotes the standard 2-norm on 
R M ). We denote this class of controls by V. The system 
equation contains a set of right invariant vector fields 
—iH\, —1H2 ■ ■ ■ , —HImi which correspond to the set of 
available one and two qubit Hamiltonians. The span 
of the set { —iH±, —1H2 ■ ■ ■ ,—iHm} (and all brackets 
thereof) is assumed to be the Lie algebra g of the group 
G 1 . Under these assumptions it follows from [TDl Prop. 



Note that we use the convention from mathematics where ele- 
ments of the Lie algebra are skew Hermitian. This is also con- 
sistent with the fact that the Hamiltonians arc Hermitian. 



3.15] that the time to move the state, from the iden- 
tity element to any other point on the group, is bounded 
(and hence, the minimum time to move between any two 
points on G is finite). Given a control signal v and an 
initial unitary Uq at time r the solution to Eq ([!]) at time 
t is denoted by U(t; v, r, Uq). 

The control problem involves generating a desired state 
Uq of the system in Eq starting from the identity 
element. By time reversal of the dynamics Eq ([I]) it can 
be seen that this is equivalent to the problem of reaching 
the identity element starting from Uq. The optimal cost 
function for this control problem is given by the geodesic 
distance 

tu (v) 

Cq{Uq)= inf { / Jv(s) T Rv(s) ds\, (2) 

VGVg 1 J > 



V s :={«(•) I HI =1, 

Vk : [0, 00) — > K is piecewise continuous 
for all k} , 

where tu Q (v) is the time to reach the identity starting 
from Uq and is defined by 

t Uo (v)=inf{t>0 : U(0) = U o ,U(t)=I, 

dynamics in ([T]) } . (3) 

This time is taken to be +00 if the terminal constraint 
U(t) = I is not attained. 

The diagonal, symmetric and positive-definite weight 
matrix R in Eq ^ reflects the relative difficulty of gen- 
erating each element of the control vector. For instance, 
on 5(7(4), the two body unitary direction <r x ® <J X may 
be weighted more than the single body unitary <j x ®I (as 
it is often harder to manipulate the former than the lat- 
ter). The symbols a x , a y , <r z denote the standard Pauli 
matrices. 

One approach to optimal control problems such as in 
Eq ([2| is the dynamic programming method developed in 
the 1950's by R. Bellman (see [H]). It relates the optimal 
cost evaluated at an initial time r at a point Uq to the 
optimal cost evaluated at a point U(r + t; v, r, Uq) that 
is reached after applying the control signal v for a time 
duration of length t. This relation takes the form 

{r+t 
[ Jv(s) T Rv{s) ds + 
J 
r 

C r+t (U(r + t;v,r,U ))\, (4) 

for all points Uq G G and is termed the dynamic pro- 
gramming equation (DPE). For any function ip this dy- 
namic programming relation can be expressed as 

C r (Uo) — SV,(r+t)(C(r+i))[^o], Uq £ G (5) 
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where 



{r+t 
/ J v(s) T Rv(s) ds + 
J 
r 



<p(U(r + t;v,r,U )) 



(6) 



Note that for ease of notation we will denote the operator 
S r .(r+t)(') as St(') when the value of r is clear from the 
context. 

In order to completely characterize the solution to the 
DPE Q we require boundary conditions given by 



C(U) 



= 0, U = I 

>0, UeG\{I}. 



(7) 



These conditions reflect the fact that if we start at the 
identity I, then the cost to reach the identity is zero. 
Furthermore, if Uq ^= I, then a non-zero amount of time 
is needed to reach the identity as the control values are 
bounded. 

The solution (i.e. the optimal cost function) to the con- 
trol problem can be obtained by solving a specific partial 
differential equation (a differential version of the DPE) , 
termed the Hamilton- Jacobi-Bellman (HJB) equation [5], 
given by 



H(U, DC) — 0, UeG\{I} 



(8) 



with 



M 



H(U,p):= sup \p- -i{J2 v kH k }U 



IMI=i 



fc=i 



with the boundary conditions in Eq Q. The DPE/HJB 
equation encodes considerable information about the 
problem, and can be solved for the optimal cost func- 
tion. The DPE can then be used to construct /verify op- 
timal control strategies via the verification theorem, [5J 
Sec. 1.5]). 

Various approaches exist to obtain the solution of the 
HJB equation ([8]). The most common are grid based 
methods [H 1121 1 1 3] which require a mesh to be gener- 
ated over the state space. Due to the COD which leads 
to infeasible memory and time requirements these meth- 
ods are unsuitable for larger dimension systems. Hence 
alternative approaches to this problem are required. 



III. THE REDUCED COMPLEXITY 
ALGORITHM 

Recently, a class of algorithms that are not subject to 
the curse of dimensionality were introduced in 17, 81 [T3] to 



solve first order HJB equations in Euclidean space. This 
method, termed the max-plus curse of dimensionality 
free approach, to solve the dynamic programming equa- 
tions involves a propagation of the solution of the HJB 
equation forward in fixed time steps without discretiza- 
tion in the spatial dimensions. The dramatic speed up 
of this approach stems from an invariant structure that 
the cost function possesses, which is preserved under the 
above propagation. This invariant form helps reduce the 
amount of information that must be stored while solving 
the control problem. In this section we introduce this 
method and describe how the invariant structure arises 
from the cost function. We then apply this method to 
obtain a numerical procedure to determine the solution 
to the control problem. 

The COD-free max-plus theory |8 to obtain and use 
an invariant form of the cost function does not currently 
deal with cost functions containing terminal constraints 
(as is the case in this problem where the trajectory must 
reach the identity element such as in Eq Q). Hence we 
formulate a relaxed version of the problem, that can be 
solved via this theory. 



A. Relaxation of the optimal cost function 

One possible relaxation of the cost function in Eq ^ 
proceeds by introducing a fixed terminal time T and a 
terminal penalty cost to yield the expression 



V a e (U ) = ferf \ I \jv{t) T Rv{t) dt + 



: <j>(U(T;v,s,U )) 



(9) 



where </>(•) is a real valued non-negative function (that is 
zero only at the identity element). This function penal- 
izes terminal states away from the identity. 

The extended control set V| above, is defined as 



V g \J{v = 0}, 



(10) 



where {v = 0} denotes the control signal that is identi- 
cally zero. We note that due to the fixed time horizon 
in Eq ^ this extension to the control set ensures that 
once the target set is reached, the cost function does not 
increase further. 

In this article we take the terminal cost to be of the 
form 



<KU) 



= tr[(/ -U)x(I- U)^} 

= tv[2J-U- U% 

= 2tr[J] -2Re(tt(U)). 



(11) 
(12) 



where 'tr' denotes the trace operation and Re(-) denotes 
the projection onto the real axis. The above relaxation 
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is valid since, as the penalty increases, the aproximation 
Vq(U) converges to V(U) for all points U G G (for a 
sufficiently large time horizon T). 

We recall that due to the assumptions outlined in 



Sec. II B the time to move between any 2 points in the 
group G is bounded. The minimum time to move from 
Uq to I is denoted by 



tu = wf e [tu («)] • 



We define 



T := sup [t Uo ] 

UoEG 



(< oo) 



(13) 



(14) 



to be the maximum value of this time over all points in 
the group. Hence choosing a T to be T in Eq ([£]) ensures 
that it is applicable for any initial point C/q G G. 

We now propagate the relaxation V* , that satisfies the 
DPE in Eq @, in time steps of r. Let T = Nt, for 
some N G N. Using the notation 

vjj(u) = vs_ kT (u) 

and Eq (JBJ we rewrite Eq ^ as the propagation of the 
function V t forward in time through a time interval r 
due to the action of the operator S T [-]: 



V£ +1 (U) = S T [V k *(-)](U), 



(15) 



where S T is the form of the operator S T which uses piece- 
wise constant controls over each time step of duration r. 

Hence we repeatedly apply S T [-] to move towards the 
desired value function Vq(-) ( = V^(-)). 

From Eq j9j, ^ it follows that for all U G G, 



v T Rv) x r+ 
V k \U{r;v,kr, CO)}, 



(16) 



which indicates the action of the DPE operator 
<§fcr.(fc+i)r[']- We now describe how the properties of this 
operator lead to the invariant structure of the optimal 
cost function Eq ([9]), and the implications resulting there- 
from. 



given 'ft' (ft G [1,2, . 
be written as 



TV — 1]), the cost function V k can 



V£(U)=c k + P k (U). 



(19) 



function becomes 



From Eqns ( 15 ) and (|16[) , after one time step r the value 



n +1 (U ) = S T [V£](U ) 
= min{(v v t j 

P k (U((k + l)r;v,kT,U ))}, 



min{(V v T Rv )r + c k + 

ueV? 



Cfe" 



= min|(v v T Rv ) t 

P k (V[(k+l)T,kT,v]-U )}, 

:= mm p v k (U Q ), 



where VC/q G G 



pl(U ) :={(Vv T Rv)r + c k + 

P k m(k + l)r, kr, v]-U ), 
*[t,s,v}-U Q :=U(t;v,s,U ). 



(20) 

(21) 

(22) 
(23) 



(24) 
(25) 



In the equations above, \l/(-) denotes the propagator for 
the system dynamics in Eq ([lj under the action of a 
control signal v. 
Hence 



T fe e +1 ([/ )= Cfc+1 TP fe+ i(f/ ), 



where 



c fc + (\ / v T Rv )r, 

Pk(*[{k + l)T,kT,v]'U ), 



Pk+i(U ) ■-- 
where v = argmin \p k (Uo) 



(26) 

(27) 
(28) 
(29) 



By the principle of induction, from Eqns. (17), (19), 



(26) it may be seen that S T [-] preserves the structure of 
the cost function. This invariance of the structure is a 
key aspect of the class of techniques introduced, as it 
helps obtain the optimal cost function at desired points 
without having to discretize along the spatial dimensions. 
This optimal cost function for any point Uq is 

VMUo) = u min min pl(U ), (30) 
ke{i,2...N} »en V T 



B. Invariant structure of the optimal cost function 



The terminal cost in Eq ( 11 ) can be written as 



where 



4>(U) = c + Pq(U), [/ G G, 



c = In, P (U) = -ti[U] - ti[U*}. 



(17) 



(18) 



Hence Co,-Po encode the initial costs corresponding to a 
control I (i.e. no control action). Assume that for a 



where Y[ denotes the ft fold product of the control set 

V%. Hence once a computationally efficient parameteri- 
zation of the control signals in terms of the set of p values 
is obtained as described above, Eq ( 30 ) easily yields the 



cost function. We note that the computation of P k and 
^[t, s, v] ■ Uq can be performed efficiently as they can be 
reduced to matrix multiplications and trace operations 
on matrices. The generation of a set of parameters, by 
using the invariant structure of the cost, and its applica- 
tion to determine the optimal cost function is the essence 
of the max-plus COD free technique. 
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From Eqns ( 23 ) , ( 24 1 it is clear that during each time towards the minimization of the cost in Eq (31). From 



step there is an increase in the size of the number of 
candidate controls to be considered during the minimiza- 
tion. Specifically, the number of elements p v k [-\ that re- 
sult from each P k is the cardinality of the control set V^. 
Thus, as in the COD free method in Euclidean space, 
due to the avoidance of spatial discretization the prob- 
lem is free of the growth in dimensionality arising from 
spatial terms; however due to the structure of the quan- 
tum control problem there is now a geometric growth in 
complexity. The details of this growth and methods to 
reduce its impact are now described. 



IV. CONTROL SPACE GROWTH AND 
PRUNING 

For the purpose of implementation let the control space 
be discretized as follows: the control signal is held con- 
stant over each particular time period of duration r. Fur- 
thermore atmost one component of the control vector is 
set to a value of 1 over any time period, while the others 
are kept at 0. This class of control signals is denoted by 
V^. Note from Eq (22) that after each time step there 



is a factor of [#(V^)] growth in the number of control 
sequences in the set to be considered, where #(A) indi- 
cates the cardinality of a set A. Hence after N time steps 
the number of possible control sequences is [^(V^)]^. 

To manage this growth we introduce a selective re- 
moval (termed pruning) of some of these control se- 
quences To describe this pruning procedure we first intro- 
duce the required notation. The set of control sequences 
of length k (i.e. k time step sequence) is denoted by 
Afc. The set of all such control sequences of all possible 
lengths {1,2 N} is 



N 



A 



fc=l 



As indicated in the previous section there is a cost 
function p\{-), associated with each control A € A such 
that the discretized cost function V^(U) = V Q (-) can be 
determined for any point U G G by 



mm p\(U), 
AeA 



(31) 



To decide upon pruning some of the control sequences 
Pa(") we first determine its contribution to the minimiza- 
tion of the cost function. A control A (and the corre- 
sponding function p^) contributes to the minimization of 
the cost function Eq ( plj ) iff 

3U G G such that P\(U) < p\(U), VA G A, (32) 

where A denotes the set of control sequences A \ {A} 
which are different from the sequence A. This idea can be 
used to measure the contribution of any control sequence 



|14) one such function that quantifies this contribution is 

C > 0, VA G A, 



W(X) = max{h((,U) \p x (U) 

U e G, Cg m}, 



where 



Px(U)=px(U)- Pl (U) VAgA, 
h(C,U) = C, (eM, 



(33) 



(34) 



If W(X) < 0, then p^ never achieves the minimum for 
any point U in G, and consequently, the control sequence 
A can be pruned without any effect on the cost function. 
For those A such that W(X) > 0, pruning would remove 
control sequences that do contribute to the minimiza- 
tion. However this is unavoidable in order to manage the 
growth in computational resources required. Therefore 
to reduce the errors in the optimal cost function aris- 
ing from this pruning, we selectively eliminate control 
sequences with relatively small values of W. This min- 
imizes the impact of pruning on the optimality of the 
resulting control strategy. 

There are several numerical methods [121 US] that can 
efficiently solve pruning problems of the form in Eq ( 33 1 . 



The more involved mathematical details of the procedure 
will be addressed in a subsequent article. By adjusting 
the upper limit on the number of control sequences stored 
in each set A&, we may arrive at an acceptable tradeoff 
between speed and accuracy. This approach has enabled 
a dramatic improvement in the time required to solve 
problems with 2-qubits (Sec. [v]), while using standard 
computing resources 



A. Description of computational complexity of the 
algorithm 

We now compare the complexity of the algorithm out- 
lined in this article with that of mesh based solution 
methods such as in [3J. In the reduced complexity 
method without pruning, the computational complexity 



grows as 



M 



K 



(35) 



where M. is the number of elements in the control set 
and K is the number of time steps in the simulation. For 
the assumptions on controllability to hold, at-most M. (G 
0{n 2 )) directions of control (i.e. control Hamiltonians) 
are required. Hence, from Eq ( 35 ) , the complexity of the 



algorithm (without pruning) for a simulation of K time 
steps is 



[0(n')] K =poly{n), 



(36) 



where poly(n) denotes a polynomial in n. With pruning, 
the complexity growth depends on the storage limits cho- 
sen in the pruning process. Hence there exists a tradeoff, 



G 



influenced by these storage limits, between the accuracy 
of the solution and the growth in complexity of the pro- 
cedure required to obtain it. 



Eq (31) for the gate synthesis problem on the two qubit 



In contrast the computational cost in Eq ( 36 ) , for mesh 



based methods with T mesh points along each dimension, 



r^- 1 ! x 0(T). 

The two terms in this expression arise from the number 
of spatial dimensions and number of iterations required 
respectively. 

An important application of the approach described 
herein is that, given a fixed time horizon T it is possible 
to efficiently (with respect to the scaling of n) check if 
a desired gate U can be synthesized within this time. 
The complexity of the algorithm to perform this check 



this would be [0(; 



,2\lK 



(without pruning) where K is 



the number of iterations in the algorithm (which is fixed 
for a given T). 

Thus it may be observed that the COD-free approxi- 
mation technique offers a potentially large order of mag- 
nitude reduction in computational complexity. 



V. EXAMPLE ON 51/ (4) 

We now proceed to apply the theory introduced, to an 
example on SU (4). The dynamics for this system is given 
by Eq with M = 5 and a control set generated from 
Hamiltonians of the form 6 {7 ® cr x , J ® <r z , a x ® 
I, a z <g> I, a x <g> a z } i.e., a set of four 1-body terms and 
one 2-body term. The associated control directions are 
sufficient to generate the entire Lie algebra su(4) , thereby 
ensuring controllability. 

To help highlight the performance improvements of the 
methods introduced herein, we note that a grid based so- 
lution approach with a conservative mesh of 50 points 
in each of the 15 dimensions of the space SU(4), with 
a total of 20 iterations over the space and an estimated 
time of 0.0001 seconds to propagate the cost function via 
value iteration [3] at each point in the mesh would require 
1.69 x 10 19 hours and a few terabytes of memory. Using 
the reduced complexity theory described in the previous 
section, this problem was solved in 15 hours to yield a 
solution for the final time horizon problem, with a hori- 
zon T of 4 seconds and a discretization step size r of 
0.2 seconds (i.e., 20 propagation steps). The simulation 
was carried out on a standard desktop computer with- 
out any exhaustive efforts to optimize the code. Hence 
there is a strong potential for further improvements to 
this procedure. 



A. Simulation results 

The simulation results obtained from the reduced 
complexity technique provide the optimal cost function 



system. In order to visualize the cost function on the 
group, we require a mapping between points on the group 
and points in Euclidean space (as the latter can be easily 
plotted via conventional graphs). For this purpose we 
make use of the exponential map from the literature on 
differential geometry jTTJ [TS] . This is an onto map that 
takes points in the Lie algebra (the tangent space at the 
identity element) to points in the group. As the Lie alge- 
bra q is isomorphic to the Euclidean space, we can thus 
obtain a function in the Euclidean co-ordinates at points 
of interest. The exponential map acts on the algebra of 
any matrix Lie group as follows 



exp(A) 



AJ 



VAeg. 



(37) 



For the connected group SU(2 n ) the exponential map of 
the algebra (exp(jj)) generates all of the group (G)[T71 
Thm 4.6], thereby ensuring that a valid visualization can 
be generated for all points in the group. 

Thus, in order to visualize the cost function on G, 
we evaluate the optimal cost function at a desired set of 
points on this group that correspond to a 2-dimensional 
slice of interest in g. The plots obtained indicate the 
approximate optimal cost function at points chosen in 
the plane of interest (on the algebra) e.g., the a x ® a x 
vs Oy ® a y plane (used in this article) . The value of the 
cost function is mapped to the shading used, in order to 
illuminate the behavior of the function. For instance in 
Fig[T] regions of darker shading indicate unitaries which 
are easier (lower cost) to generate while the lighter areas 
show gates which are costly to synthesize. 

In order to understand the effects, on the cost func- 
tion, of the difficulty in synthesizing the 2-body unitary 
compared to the 1-body terms (which we denote by the 
ratio r of the cost of one body interactions to that of the 
available two body interaction a x (g> <r x ) , simulations were 
performed that varied r (using the R matrix in Eq ([2])) 
starting from slightly less than one and proceeding upto 
a value of one-third. 

The generation of unitaries in the a y ® o~ v direction 
(which is not directly accessible) requires the alternating 
application (i.e., bracketing operation) of o~ x (8> & x and 
one of the available single body Hamiltonians. Hence if 
r = 1/1.3, then the cost of moving along the <r y <8> cr y 
direction is substantially larger than that along the 
one and two body control directions, leading to an 
elliptical shape of the level set (as shown in Fig. |2j, 
where the cost increases faster in the o~ y ® o~ y direction 
compared to the o~ x ® a x direction. However as the 
cost imposed on moving along the <r x ® <r x direction 
is increased, the relative cost of moving along a y <£> a y 
compared to o x ® a x decreases, thereby leading to level 
sets that become more circular (Fig. This agrees 
with analytical results such as in [19]. Note that there 
are two axis of symmetry in these plots, namely both 
o~ x ® & x and <j y ® o~ y . This arises due to the symmetry 
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FIG. 1: Plot of the cost function in the plane a x ® <? x vs 
a y <g> cr v . Here the cost of applying any available 1-body 
Hamiltonian is taken to be only slightly less than the cost of 
generating the 2 body Hamiltonian a x ® a x (i.e., r = 1/1.3). 
The lighter areas indicate unitary operations that are harder 
to synthesize. 




* -jr/2 1/2 

a (9 a 

x x 



FIG. 2: Top view of the plot in figure [T] 

in the cost function equations |2]) and Q with respect 
to the application of controls along cither the positive or 
negative direction of the available control Hamiltonians 
in the system equation 



VI. CONCLUSIONS 

In this article we have demonstrated a reduced com- 
plexity method may be used to obtain substantial im- 
provements in computational speed in solving optimal 
control problems arising in closed quantum systems. This 



FIG. 3: Plot of the cost function in the plane a x ® a x vs 
o y ® o v with r set to 1/3. 
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a 8 c 

FIG. 4: Top view of the plot in figure [3] 

technique was used to obtain a numerical solution for an 
example gate complexity problem in a 2 qubit system. In- 
stead of the curse of dimensionality in the spatial dimen- 
sion, we now have a much more manageable growth in 
dimensionality that depends on the number of elements 
in the discretization of the control set. 

The approach outlined in this article can deal with a 
very general class of problems and can, in principle, be 
used for systems with any number of spins. Furthermore, 
the methods described can be extended to other systems 
of interest. For instance, a control problem on a sys- 
tem with drift can be solved under the current system 
framework (Eq ([I])) by taking the cost of moving along 
the negative direction of the drift term to be much larger 
than that in the positive direction. 

At present, the techniques outlined yield preliminary 
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solutions whose error bounds, rate of convergence and 
other properties must be determined via further research. 
It is hoped that the method introduced (and refinements 
thereof) would enable the accurate solution of problems 
on spin systems of larger dimensions than has been pos- 
sible until now. 
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